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LIST  OP  SYMBOLS 


SUPERSCRIPTS 


n  refers  to  the  time  index, 

o  refers  to  time  zero. 


SUBSCRIPTS 


k  refers  to  the  R  space  index. 

I  refers  to  the  Z  space  index, 

a  refers  to  charge  a. 

b  refers  to  charge  b. 


GREEK 


a  =  an  adjustable  constant  to  steepen  shock  fronts, 
p  =  density. 

ROMAN 


c  =  sound  velocity. 

E  =  specific  internal  energy. 

Sr  Sr 

55  3z 

J  = 

Sz  Sz 

3R  Sz 


Sp  Sp 

55  5z 

k  = 

Sz  Sz 

55 


SP  SP 

55  5z 

l  = 

Sr  Sr 

35  5z 

p  *  pressure. 

P  =  p  +  q. 
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q  ^ 


r  = 


R  = 

m  = 

t  = 
At  ■ 

u  = 


v  = 


w 


o 


z  = 


AZ  = 


la °f  i 

p 


if 


<  0  • 


if 


if 


>  °  , 


Eulerian  coordinate  perpendicular  to  the  z-axis  of  symmetry 
initial  radius  of  charge  a. 


initial  radius  of  charge  b. 

Lagrange  coordinate  perpendicular  to  the  Z-axis  of  symmetry 

increment  of  Lagrange  coordinate  R. 

time. 

increment  of  time. 

dr(tfRjZ)  =  r  component  of  velocity. 


=  z  component  of  velocity. 


St 


o  o  To 
p  r  J  . 

Eulerian  coordinate  of  cylindrical  symmetry. 


initial  axial  distance  of  center  of  charge  a. 


initial  axial  distance  of  center  of  charge  b. 
Lagrange  coordinate  of  cylindrical  symmetry, 
increment  of  Lagrange  coordinate  Z. 
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INTRODUCTION 


Axi- symmetric  unsteady  compressible  flow  problems  arise  extensively  in 
military  technology.  Spherical  blast  waves  and  their  interaction,  explosions 
in  vertically  variable  atmosphere,  flow  past  axi-symmetric  bodies,  high-speed 
flow  in  axi- symmetric  hypervelocity  guns,  shaped-charge  gas-metal  jet  flow-- 
all  are  of  importance.  The  exact  solution  of  the  equations  governing  these 
flows  which  satisfies  the  appropriate  initial  and  boundary  conditions  requires 
solving  a  system  of  non-linear  hyperbolic  partial  differential  equations  in  a 
two-space  and  one-time  coordinate  system,  i.e.,  in  three  generalized  dimensions. 
Presently,  such  exact  solutions  are  not  obtainable.  The  method  of  character¬ 
istics  by  which  existence  of  solutions  is  shown  opens  the  way  for  a  numerical 
procedure,  but  this  method  requires  a  series  of  machine  programs  that  esca¬ 
lates  in  complexity  as  the  shock  and  characteristic  discontinuity  surfaces  in 
the  flow  field  interact  and  multiply  the  number  of  continuous  zones  bounded 
by  discontinuity  surfaces. 

Von  Neumann  and  Richtmyer  proposed  a  method  of  numerical  integration  that 
changes  the  flow  equations  to  the  parabolic  type.  Although  the  solution  is 
only  approximate,  the  overall  or  essential  nature  of  the  solution  is  preserved; 
i.e.,  curves  with  momentum  and  energy  discontinuities  are  warped  into  smooth 
but  rapidly  changing  curves.  (Stream  surface  material  discontinuities,  however, 
are  not  altered.)  Therefore,  the  numerical  procedure  does  not  require  the 
separate  treatment  of  every  continuous  zone,  and  numerical  solution  with  today's 
computing  machines  is  quite  feasible. 

This  is  the  method  to  be  examined  for  calculation  of  the  interaction  of 

the  blast  waves  from  two  spherical  charges  of  initial  radii  r°  and  r°,  since 

a  d 

in  this  interaction  the  flow  field  is  replete  with  discontinuity  surfaces  that 
complicate  the  problem. 
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Shock 


FIGURE  1.  ELAST  INTERACTION 

While  this  work  is  the  theoretical  extension  of  the  experimental  work  done 
in  References  1,  2,  and  3  on  project  CLAW,  with  small  modifications,  the 
analysis  can  be  extended  to  the  other  military  problems  previously  listed. 
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PLOW  EQUATIONS 


The  continuity  equation  in  particle  or  Lagrange  cylindrical  coordinates 
for  axi- symmetric  flow  is 


(1) 


(prJ)  =  o,  J  = 


Sr 
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or 


(2) 


T  o  o  To  o 
prJ  =  p  r  J  =  w  . 


The  momentum  conservation  equations  are 
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where 


f  Sr 

St 


(4)  1 


Sz 

St 


The  energy  conservation  equation  is 


/r\  SE  P  Bp 

(5)  St  •  “2  St 


li 


where,  in  accordance  with  Reference  h,  P  is  defined  by 


(6)  P  -  p  +  q  , 


and  q  represents  some  artificial  dissipative  term.  In  accordance  with 
References  5  and  7  we  put 


(7) 


q  = 


>  o  , 


The  conservation  equations  above  are  supplemented  by  the  thermodynamic 
equations  of  state 

(8)  P  =  p(E,p)  , 

c  =  c(E,p)  , 

which  are  assumed  to  be  known. 

The  conservation  equations  above  are  hydrodynamic ally  exact  for  q  =  0; 
but  for  q  ^  0,  the  momentum  and  energy  equations  are  modified,  while  the 
mass  equation  remains  intact. 

Equations  (2)  thru  (8)  are  ten  equations  with  the  ten  unknowns 
(r,  z,  u,  v,  E,  P,  p,  q,  p,  c),  which  are  to  be  determined  as  functions  of 
(t,  R,  Z). 

INITIAL  AND  BOUNDARY  CONDITIONS 

o  o  o  o  o 

Two  explosive  charges  of  radii  r  and  r.  with  centers  at  z  and  z.  ■  -z 

a  d  a  d  a 

along  the  Z-axis  (Figures  2  and  3),  are  initiated  simultaneously  at  their 
centers.  Assuming  equilibrium  flow  and  steady  state  Chapman- Jouget  propagation 
of  the  detonation  front,  the  initial  explosion  state  can  be  calculated  either 
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by  the  Taylor  progressive- wave  approximation  (Ref.  8),  or  by  the  numerical 
scheme  to  be  subsequently  described.  In  either  case,  the  equation  of  state 
of  the  explosion  products  in  chemical  and  thermodynamic  equilibrium  must  be 
specified. 

If  the  detonation  process  is  to  be  calculated  by  the  method  outlined  in 
this  article,  considering  that  the  media  ahead  and  behind  the  detonation 
front  are  not  the  same,  the  wave  fronts  must  be  treated  as  boundaries  moving 
with  the  constant  Chapman- Jouget  velocity.  The  wave  fronts  form  cones  in 
the  R,  Z,  t  space  (see  Figure  2).  Behind  the  detonation  fronts  the  calcu¬ 
lations  use  the  equation  of  state  of  the  explosion  products.  Ahead  of  the 
fronts  are  the  undetonated  portions  of  the  charges. 

When  the  detonation  process  is  complete  and  the  two  explosion  gases  expand 
into  the  surrounding  ambient  air,  shock  waves  are  formed  and  propagate  in  an 
outward  direction.  These  shocks  eventually  collide,  and  the  two  flows  interact. 
The  primary  area  of  interest  is  in  the  nature  of  the  interaction. 

On  each  t  =  constant  plane  (Figure  4)  after  the  shock  formation,  two  types 
of  domains  must  be  considered  --  the  explosion  gas  and  air  domains.  Each  has 
its  own  equations  of  state  (8).  Although  shock  surfaces  are  moving  boundaries 
of  discontinuities,  the  momentum  and  energy  equations  modified  by  the  higher 
derivative  dissipative  term  no  longer  permit  this  type  of  discontinuity  sur¬ 
faces  to  exist.  In  essence,  such  boundaries  are  smoothed  out  and  no  longer 
distinct. 

In  the  Lagrange  space  (Figures  2  and  4), each  explosion  gas  boundary 
remains  a  fixed  circle  in  each  t  =  constant  plane.  The  equations  of  these 
circles  are 


(9) 


R2  +  (Z  -  z  f  =  (r  )2  , 
\  ay  v  a7  * 


VR2  +  (Z  -  *b)2  =  (rb)2 


Across  these  circles  we  assume  dynamic  equilibrium  and  laminar  flow  so  that 
pressures  and  velocities  normal  to  the  stream  surfaces  are  continuous.  The 
stream  surfaces  remain  as  characteristics  across  which  thermodynamic  variables 
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FIGURE  4.  t  =  CONSTANT  PLANE 
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may  be  discontinuous,  because  the  mass  equation  is  unaltered  in  the  modified 
system  of  conservation  equations.  Thus,  the  chemical  composition  and  equations 
of  state  change  abruptly  across  the  explosion  gas-air  boundary. 

COMPUTATIONAL  SCHEME 


The  computational  scheme  for  points  interior  to  boundaries,  i.e.,  points 
other  than  explosion  gas-air  boundaries,  will  be  that  adapted  to  this  particular 
problem  from  Reference  5.  The  stability  of  this  scheme  is  shown  in  Reierence  6. 


Let  k,  /,  n  be  indices  corresponding  to  the  coordinates  R,  Z,  t  (Figure  5), 
Assume  that  on  the  t=  0  or  n  =  O'th  plane  all  quantities,  including 


r  \  1  JL  =  p  _1  _1  r°  1  _1  J°_l  _1 


1  /  o  *  ,  o  ,  ,  .  o  ,  t  o  . 

(v,  X  1  +  W.l  1  +  V  1  1  +  w,  1  1  )  , 


(10) 


< 


w,  .  =  r  J-  -  1  +  w.  1  .  1  +  v.  1  M  1  +  w,  1  .  l 

k,i  4  k-^,  ^"2  ^~29  ^"2  k+2>  ^*2 


^,1  -  *  <^4,  4 +  «k4.  i4 +  <4 +  **4.  «4  >  ’ 


2,  *.2 


2'  2 


"S’  inr2 


AZ, 


k,i  *  ¥  <AZ*4,  <4  *  ffik4,  ,4 +  fflk4,  <4  *  ^4,  «4  )  - 


2*  2 


2'  2 


r2»  4T2 


are  known  at  the  pre- determined  points 


f  Z 


and  R,  „  •  Z,  . 
k,  i 7  k,i 


(This  grid  spacing  in  the  R,  Z  space  may  require  changing  at  some  time  tn, 
depending  on  the  gradients  of  the  solution.)  Also,  assume  that  on  the 
successive  (n-^)th  and  (n-l)th  constant  time  planes,  the  following  quantities 
have  been  computed  for  all  values  of  the  indices  k,  i  up  to  k,  i: 
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(n-|)th  plane 

(At)n4  , 

n  1  „  1 

U.  1  1  *  V.  1  .1 

K-j>  i-15  K-^,  *-75 

(n-l)th  plane 

v1  , 


n-1 

Pk,i 


z 


n-1 

k,i 


y 


n-1  n-1  n-1 

J  1  l  9  p  l  i  9  P  i 

k-5,  k-^,  k-^, 


n-1 


n-1 

C,  1  1  * 
k-^,  i-jy 


The  computational  network  sequence  is  chosen  such  that  on  each  constant  n 
plane  k  varies  along  constant  i  strips. 


On  the  nth  plane,  r  and  z  are  approximated  using  the  velocity  definitions  of 
equation  (4). 


(n) 


r 

5 


n 

rk,  l 


n-1  n  -  | 
rk,i  +Uk,i 


1 

2 


n 


z 


k,i 


n-1  n  -  ^ 
*  vt,i 


n 

At 


1 

5 
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At  half  points  r,  z  are  obtained  by  averaging  the  values  at  the  neighboring 
whole  points. 


n  1/  n  ,  n  n  ^  n  x 

*  1  ,  Is  U(rk-1,/-1  +  rk,/-l  +  rk-l,£  +  rk,f} 

k  "  -  2 

n  1/  n  n  ,  n  .  n  x 

z  1  1  =  Trzk-l,i-l  +  zk,/-l  +  zk-l,i  +  zk,i' 

“  3 9  *  **  5 


From  these  values  the  gradients  of  r,z  are  approximated  by  the  difference 
quotients 


which  give  for  the  Jacobian 
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The  second  and  third  equations  in  (l8)  will  differ  for  air  and  explosion  gas. 
P  is  approximated  by 
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roll-1 


With  the  quantities  computed  previously, 


the  following  can  be  calculated 
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Minimizing,  the  following  is  obtained: 


If  k  -  1,  l  -  1, 


D  =  MIN 
1  1 
2  '  2 


If  k  =  1,  l  >  1, 


(  Bn  ,  Cn  )  . 

1  1  11/ 

\  5*2  2 ’ 2/ 


Dn  =  MIN  /  Bn 


1  *  1 
"  2 


/  B  C“  D  v 

(1,1  1*1  max(k*-  i  ),i  -  li  )  . 

V  2ji  *  2  »  2»i  ■  2  »  2  2  / 


If  k  >  1  , 


d;  ^  *  MIN 

k  -  ?>>l  -  ^ 


/  Bn  .  .  ,  Cn  -  ,  Dn  \ 


(Note:  The  calculation  above  assumes  that  k  increases  on  constant  l  strips.) 


Time  increments  are  calculated  from 


n+|  /  n-in  \ 

t  =  MIN  (  1.4  At  ,  D  ) 

V  max(k  -  ^),  max(  t  -  i>) J 

,n  1  /  A+n  ‘  2  “  +  2  ) 

t  =  J  (  At  +  At  /  , 


and  the  time  from 


n  -  1 


+  At 


1 
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n  +  *  n  -  2  n 

t  *  t  +  At 
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k- 1,2-1  “k-1,1-1 


Prom  these  results,  K  and  L  are  evaluated  as: 
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k-1,1-1 


(26) 
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■8Pvn 
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<s> 

SPv11 
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The  momentum  conservation  equations  (3)  give  the  accelerations 


,du* 

(St) 
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from  which  the  velocities  are  derived: 
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This  completes  the  evaluation  of  all  points  away  from  the  explosion  gas-air 
boundaries. 


Each  network  point  k,  £  is  examined  to  determine  if  it  crosses  any 

boundary.  Proceeding  on  successive  constant  Z  or  i  lines  from  R  =  0  to 

R  >  0  (Figure  4),  starting  from  a  negative  Z  <  (z°  -  r°)  to  some  positive 

Z  >  (z°  +  r°),  the  R  =  0  point  has  entered  the  domain  of  the  explosion  gas  b 
a  a 

whenever 


rb  i 


0,1 


.  o  .  O 

i  strb 
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and  has  entered  the  domain  of  the  explosion  gas  a  whenever 


o 

z 
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<  Zo,£ 
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o 

z 

a 


+  r 


o 

a 


If  the  R  =  0  point  is  inside  the  explosion  gas  b,  along  each  Z  =  Z 

K  ,  jb  O  y  l 


=  constant  line,  the  inequality 


A 


,i  - ^ 


°\2  /' 7  °t2 

b>  '  <Zo ,t  -  zb> 


J  <  0  , 


obtained  using  (9)  must  hold  as  R^  increases.  If  the  inequality  holds,  R 

is  inside  the  explosion  gas  b;  otherwise,  R  is  in  air.  Similarly,  if  the 

R  =  0  point  is  inside  the  explosion  gas  a,  along  each  Z  =  Z  =  constant 

K ,  £  0 ,  £ 

line,  if 


!  ,  -  (r°)2  -  (Z  -  z°f  <  0 

i,£  '  a  v  o,£  a'  - 


is  satisfied,  R  is  inside  a,  and  outside,  if  the  inequality  fails. 


At  any  stage  where  R  traverses  an  explosion  gas  boundary  into  air,  J,  K, 

L  are  evaluated  by  taking  finite  difference  quotients  across  this  discontinuity 
boundary.  While  this  procedure  is  not  exact,  the  error  committed  attenuates 
with  propagation  into  the  flow  field. 
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CONCLUSION 


The  new  BRLESC  computing  muchine  of  BRL,  with  addition  times  in  micro¬ 
seconds,  can  compute  very  rapidly  the  three  dimensional  network  required  in 
the  numerical  scheme  of  this  problem.  While  the  extensive  information-hunting 
and  reading  with  the  magnetic  tape  will  slow  the  overall  computation  considerably, 
a  total  of  52,000  fast  core  memory  will  be  available  soon,  making  the  problem 
entirely  feasible  and  the  computing  time  very  fast. 


C, 


R.  C.  MAKINO 
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